#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double
generate_gaussian_noise(double variance)
{
	static double rand1, rand2;
	static bool hasSpare = false;
	static double TWO_PI = (2 * M_PI);

	if(hasSpare)
	{
		hasSpare = false;
		return sqrt(variance * rand1) * sin(rand2);
	}

	rand1 = ((double) rand()) / ((double) RAND_MAX);

	if(rand1 < 1e-100)
		rand1 = 1e-100;

	rand1 = -2 * log(rand1);
	rand2 = ((double) rand() / ((double) RAND_MAX)) * TWO_PI;

	hasSpare = true;
	return sqrt(variance * rand1) * cos(rand2);
}


double
polar_box_muller()
{
	double x1, x2, w, y1, y2;

	 do
	 {
		 x1 = 2.0 * ((double) rand() / ((double) RAND_MAX)) - 1.0;
		 x2 = 2.0 * ((double) rand() / ((double) RAND_MAX)) - 1.0;
		 w = x1 * x1 + x2 * x2;
	 } while (w >= 1.0);

	 w = sqrt((-2.0 * log(w)) / w);
	 y1 = x1 * w;

	 return y1;
}

int
main(void)
{
	double variance = 0.1;
	srand(time(NULL));

	for (int i = 0; i < 1000; i++)
		// printf("%lf %lf\n", generate_gaussian_noise(0.1), generate_gaussian_noise(0.1));
		printf("%lf %lf\n", polar_box_muller() * variance, polar_box_muller() * variance);

	return 0;
}
